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COOLING EFFECT OF THE RICHTMYER-MESHKOV INSTABILITY 

F. MohsenF, M. Mendoza^, S. Succi^ and H. J. Herrmann^’^ 


Abstract. We provide numerical evidence that the Richtmyer-Meshkov (RM) instability contributes 
to the cooling of a relativistic fluid. Due to the presence of jet particles traveling throughout the 
medium, shock waves are generated in the form of Mach cones. The interaction of multiple shock 
waves can trigger the RM instability, and we have found that this process leads to a down-cooling of 
the relativistic fluid. To confirm the cooling effect of the instability, shock tube Richtmyer-Meshkov 
instability simulations are performed. Additionally, in order to provide an experimental observable of 
the RM instability resulting from the Mach cone interaction, we measure the two particle correlation 
function and highlight the effects of the interaction. The simulations have been performed with an 
improved version of the relativistic lattice Boltzmann model, including general equations of state and 
external forces. 


1. Introduction 

Particles traveling through a compressible fluid generate waves moving at the speed of sound. Moreover, if 
the particles travel faster than the speed of sound of the medium, the disturbances in the fluid are confined to 
the so-called Mach cone. This phenomenon is very common in many natural systems, including astrophysics and 
high energy physics CHS], where relativistic fluid effects are important. The existence of relativistic shock-waves 
in the presence of density variations, leads to the appearance of the RM instability, one of the fundamental 
fluid instabilities, which occurs whenever a shock wave passes through an interface between regions at different 
densities. This instability was theoretically predicted by Richtmyer [6] and experimentally detected by Meshkov 
[7], in the non-relativistic context. The study of the RM instability is of major importance in several fields, 
ranging from high energy physics mg to astrophysics [8] and plasma physics [9] . Density variations can appear 
in the relativistic fluid whenever particles travel through the medium, due to the sweeping effect of the shock 
waves m, as well as due to external mechanisms. 

In this work, we show numerically that the RM instability may reduce the average temperature of the 
relativistic fluid. In particular, we investigate the interaction of two relativistic Mach shocks, and show that 
the RM instability arises due to this interaction (see Fig. [1]). Furthermore, we find that the appearance of 
this hydrodynamic instability leads to a decrease in the average temperature of the medium. To justify this 
finding and to single out the effect of the instability, shock tube RM instability simulations are carried out. The 
effect of initial domain temperature and density ratio on the cooling effect of the instability is also investigated. 
Since the growth rate of the instability depends explicitly on its form, the study of the instability can provide 
information on the equation of state (EoS), the same way shock waves can offer insights on the EoS 
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Thus, gaining information about the RM instability may provide a new means of studying the thermodynamic 
properties of relativistic fluids. Finally, we propose a way to detect the interaction between Mach cones from an 
experimental observable, namely the two-point correlation function. For the numerical simulations, the recently 
developed relativistic Lattice Boltzmann (LB) model [13] is extended to deal with the ideal gas equation of 
state and external forces. 

This paper is organised as follows. In Sec. [21 we explain the numerical model and the extensions needed 
to describe a relativistic fluid with an arbitrary equation of state. The results of our numerical simulations of 
the RM instability are presented in Sec. |3l and finally, in Sec. H] we discuss our results and comment on future 
works. 


2. Numerical model 

Let us start the description of the numerical model by presenting the conservation equations for relativistic 
fluid dynamics, namely 

= 0 , daN^ = 0 , ( 1 ) 

where the energy-momentum tensor is defined as = {e+ p)U^U^/c^ and the current density 

as = nU^ [14]. Here, n is the number density, p the hydrostatic pressure, e the energy density, c the speed 
of light, the shear-stress tensor, and the Minkowski metric tensor with the signature (-h, —, —, —The 
macroscopic four-velocity is (U^) = (c, u)j{u), u being the three-dimensional velocity and ^{u) = l/^/l — v?/ 
the Lorentz’s factor. The Einstein summation convention and natural units i.e., c = ks = h = are assumed 
here and throughout this paper. 

All our numerical simulations are performed using the extended version of the relativistic LB model recently 
proposed in Ref. m- This method is a numerical approach based on a minimal lattice version of the relativistic 
Boltzmann equation, which can be solved to find the probability distribution function in phase space p^llTK] . The 
equilibrium distribution for the relativistic Boltzmann equation, in the single-relaxation time approximation 
for the collision operator, is the Maxwell-Jiittner distribution function [16]. An extension to simulate high 
velocities was also proposed in Ref. [13], based on the ultra-relativistic equation of state, i.e., e = 3p. However, 
for this case, the equation for the conservation of energy and momentum is not affected by the density field m, 
and the two conservation equations become decoupled. This effect suppresses the RM instability, where both 
equations must be coupled. Therefore, we are interested in a more general ideal gas equation of state, of the 
form [TT] : 

P= (r - l)(e-n), (2) 

where F = Cp/c^, with Cp and Cy being the specific heats at constant pressure and volume, respectively. For low 
temperatures, i.e. mc^ /ksT ^ 1 , F = 5/3, while for high temperatures, i.e. me? jksT ^1 , F = 4/3. In the 
ultra-relativistic limit, by replacing F = 4/3 and considering the condition n ^ e, the ultra-relativistic equation 
of state is recovered. 

In the relativistic lattice Boltzmann methods based on the model of Marie [18] for the collision operator, 
the macroscopic variables can be calculated by solving a system of equations, corresponding to the moments 
of the equilibrium distribution [mus]. However, for the case of the ideal gas equation of state, this system 
of equations cannot be solved due to the fact that the first and the second order moments are coupled. This 
problem can be solved by using the model of Anderson-Witting m for the collision operator. According to the 
model of Anderson-Witting, the relativistic Boltzmann equation has the following form 

TT a 

P^daf =- n), ( 3 ) 

T 

where (p^) = ( 7 ,^ is the four-momentum, / is the probability distribution function, and r the single relax¬ 
ation time. The model of Anderson-Witting is based on the Landau-Lifshitz decomposition m- Hence, the 
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macroscopic variables can be calculated using the following relation m 


UaT^^ = eU^. 


(4) 


In this case, e and are the largest eigenvalue and corresponding eigenvector of T^, respectively. These values 
can be calculated numerically using the power method. The density, and the pressure would be evaluated 
subsequently, using the first order moment relation and the equation of state, respectively. 

Additionally, in order to equip the model with an ideal gas equation of states, the discretised distribution 
function should be also modified. Hence, the following term should be added to the original distribution function 
proposed in Ref. 

3(r - l)(e - n) - e 52^(3(r - l)(e - n) - e) 

(r-l)(e-n) + e 33 x 72((r - l)(e - n) + e) ’ 

where (5^0 is the Kronecker delta function. Note that, as expected, in the ultra-relativistic limit, this term goes 
to zero. For other details of the numerical model and its validation, one can consult the original article [13]. 

In the presence of disturbance traveling in a relativistic fluid, the energy-momentum conservation can be 
expressed as = 5'^, where the source term is the energy deposited by the disturbance and can be 

written in the form 1211122]: 


= 


1 


■ exp 


{- 


[X -Xjetf 


dE ^ 


(6) 


I 2<t2 

where the momentum deposition is ignored and the disturbance is assumed to travel with the velocity of light. 
Here Xjet is the location of the disturbance, where a = 0.04 and dE/dx = 7.5 are considered. Numerical units 
are used here and throughout the paper. 

Additionally, to include the external force, *5^, into the relativistic LB model, we need to calculate the 
discretised forcing term which will be added to the discretized Boltzmann equation. For the relativistic LB 
scheme in Ref. [13], and assuming the external force as {S^) = (5'^, S'), the discretized forcing term becomes: 


S, = 


2u 


, 2^3 


5t 




cts^ 

Coe p 


(7) 


where / is defined in Eq.([5]), and Wi are the discretized lattice vectors and weight functions, respectively, and 
u, Co and Ct are lattice constants which can be found in Ref. m- 


3. Results 

For the numerical simulation of the interaction between two Mach cones, two domains with 500 x 500 and 
250 X 250 cells are considered. All boundaries are taken as free outlets and the ideal gas equation of state is 
assumed with F = 4/3. The initial position of the disturbance moving along the x direction is (I//3,1//2), while 
the disturbance moving along the y direction starts at L being the length of the domain. We take 

St/Sx = 0.15 and the controlling parameter of the bulk viscosity is a = 0.25 (see Ref. [I3]). For more details on 
the numerical technique, see Ref. [23] . 

In Fig.lU we show that, after the interaction of two Mach cones, when the shock front of a Mach cone passes 
through the density variation which is caused by the other Mach cone, the RM instability starts to grow in the 
direction of the advancing shock front. In this figure, we also show that downstream the moving disturbance 
and due to the sweeping effect of the shock wave, the density decreases locally. It is worth mentioning that 
the current simulation shows only one particular case of the interaction of the Mach cones. Nevertheless, what 
causes the instability is the interaction of one Mach cone with the density fluctuation due to the passage of 
another Mach cone at earlier times. 

In order to inspect whether the current grid resolution is sufficient to capture the shape of the interaction, 
additional simulations, with the same parameters as before, but with 700 x 700 cells are reported in Fig. [2l In 
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Figure 1. Snapshots of the density profiles for two interacting Mach cones in a relativistic 
fluid at times (a) t = 180 (b) t = 360 (c) t = 540 (d) t = 720 (e) t = 960, and (f) t = 1200. 
Here Tmed = 0.3, the red and blue colours denote high and low values of the density, respectively. 



(a) (b) 

Figure 2. Snapshots of the density profiles for two interacting Mach cones in a relativistic fluid 
at time t = 1200 using lattice sizes of (a) 500 x 500 and (b) 700 x 700 cells. Here Tmed = 0.3, 
the red and blue colours denote high and low values of the density, respectively. 











ESAIM: PROCEEDINGS AND SURVEYS 


5 


H 

<1 


C/2 


0.03 


0.02 


0.01 


0‘ 


A A 


A T„,,=0.4 

□ T„,,=0.8 


oooogoooooo^ 

"a ° 

A A^ 


□ □ 


O □ □ □ □ □ 

^ ° ° 


1000 1500 2000 

Time 


Figure 3. Decrease in the average temperature of a relativistic fluid due to the RM instability 
developed as a result of the interaction of two Mach cones, for different initial temperatures, 
Tmed- Here ATins = Tq — with T^g^ being the average temperature and To the average 

temperature right after the disturbances have left the simulation zone. Thus, ATins denotes 
the decrease in temperature of the media compared with the time when the disturbances and 
shock waves had left the domain. Here the domain of simulation is 250 x 250. 


the same figure, we also show a late stage of the interaction. Despite a slight difference between the results at 
different resolutions, one can appreciate that the shape of the instability due to the interaction is well captured 
with the current 500 x 500 lattice resolution. 

Regarding the thermal behaviour of the fluid during this phenomenon, as expected, the passing disturbance 
in the medium increases the average temperature, since it deposits energy to the fluid according to Eq. 
Here, we are interested in the effects of the aforementioned RM instability on the average temperature. Thus, 
we compute the average temperature of the medium and compare it to the average temperature when the 
disturbances and shock waves have left the domain completely. In the absence of the instability, the temperature 
should remain constant because of the steady state condition. However, we see that, due to the presence of the 
instability, the temperature starts to decrease, see Fig.[3l The simulation is performed for different initial domain 
temperatures and we also observe that the decrease in the initial temperature enhances the cooling effect. This 
can be explained by realizing that relativistic effects, which are more dominant at higher temperatures, weaken 
the RM instability, as recently shown in Ref. [23]. More precisely, the linear growth rate of the instability (vf) 
is found to take the following form [23] : 


{n2-ni)khoAu 

= 7(2p + « + <i)' 

where, k is the initial wave number, Hq is the initial amplitude of the perturbation, Au and p are the velocity 
and pressure at the interface, respectively, and 62 ( 722 ) and ei (ni) are the energy (density) at both sides of 
the interface. This expression suggests that the growth rate of the relativistic RM instability decreases as the 
temperature increases, i.e. the denominator increases in Eq.([H|). This is in agreement with the results in FiglS] 
where we show that a higher initial temperature leads to a smaller decrease in the average temperature due to 
the RM instability, since the instability is weaker at higher temperatures. 

At this point, we have shown that the RM instability can appear during the interaction of Mach cones in 
a relativistic fluid, and that it contributes to the cooling of the medium. In order to single out the instability 




6 


ESAIM: PROCEEDINGS AND SURVEYS 



Figure 4. Snapshots of the density (top) and temperature (bottom) fields at t = 1260 in the 
shock tube RM instability with perturbed interface. Here, we consider the high density ratio 
^l/^m = 28 and Mar = 2.4. Blue and red colors denote low and high values, respectively. 


and to confirm its cooling effect, we perform simulations in a simplified configuration, namely the shock tube 
RM instability. The reason for choosing a simple shock tube geometry is twofold; first, this is a standard 
geometry to study the RM instability, second, it is straightforward to compare the cases with and without the 
instability. Simulations are performed on a square lattice with 1200 x 200 cells. For all simulations considered 
here, a shock wave with the velocity P = \u\/c = 0.94, travelling from right to left, is passing through a 
sinusoidal perturbation in the density, located at Xp = 1000 cells. The sinusoidal perturbation takes the form, 
Xi = Xp asin(7r/2 + ^'KyjX)^ where A is the width of the domain and a = 32 (for further technical details 
about the numerical simulation, see [23]). Note that subscripts R, M , and L refer to the right hand side of the 
shock, the region between the shock and the initial perturbation, and the left hand side of the perturbation, 
respectively. Periodic boundary conditions are considered for the top and bottom boundaries of the domain, 
while inlet and outlet boundary conditions are applied to the right and left boundaries, respectively. Other 
parameters are chosen as follows, StjSx = 0.15, a = 0.25 and F = 5/3. For comparison, the cases without the 
instability are also simulated by simply setting a = 0 (unperturbed interface), while other parameters are the 
same as the ones for the case with the instability. The snapshots of the density and temperature profiles, for the 
perturbed case with the density ratio Rl/^m = 28 and relativistic Mach number Mar = 2.4 at a late time of the 
instability, are presented in Fig. |4] {Mar = Us^{us)/ca^{cs) with Ug being the shock velocity and Cg the sound 
speed). By measuring the average temperature of the plasma on both cases, with (perturbed interface) and 
without (unperturbed interface) instability, the decrease in temperature due to the instability can be computed. 
In Fig. [SI one can notice that in the case of a perturbed interface, the average temperature is lower, where the 
cooling effect of the instability increases in time. These numerical experiments are in agreement with the results 
shown in Fig. O The results for different density ratios in Fig. [5] show that, as the density ratio increases, the 
decrease in the average temperature is enhanced. This is because at higher density ratios the instability grows 
faster, than predicted by Eq.(|8|). 

Back to the our original problem of Mach cone interaction, we explore the possible experimental observable 
consequences of the existence of this kind of interactions. To this purpose, we suggest to study the two-particle 
correlation (TPC) function [T2l|22]. Our hydro dynamical calculations provide macroscopic quantities such as 
temperature and velocity fields. Thus, to compare the hydro dynamical results with experimentally measured 
observables, a description of the conversion of the fluid into particles is needed. This can be achieved by the 
Cooper-Frye freeze-out approach [24], where the particle emission pattern is given by: 


where Pt = yvi + Py is the transversal component of the momentum of a particle, (j) is the azimuthal angle, j 
is the pre-factor for the Maxwell-Jiittner distribution function, y = 1/2 In is the rapidity and is the 
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Figure 5. Decrease of the average temperature due to the RM iustability iu a shock tube 
uumerical experimeut, for differeut deusity ratios. Here < Tunp > aud < T^er > deuote the 
average temperature wheu the RM iustability is (perturbed iuterface) aud is uot (uuperturbed 
iuterface) preseut, respectively. We have set Mar = 2.4. 


iutegral surface iu space-time. Usiug isochrouous freeze-out hypersurface, = (l,0)d^:r, we have [22] 

dN 

PTdpTd(f)dy 

j 


^ J d^xrriT cosh ^ X exp{ — ^ [ruT cosh ^ 


( 27 rj 

—prUx cos (j) — prUy siu (j )]} 


( 10 ) 


where mr = E/ cosh^. 

Upou defiuiug: 

with /o the Maxwell Jiittuer distributiou fuuctiou at T = Tmed aud r = 0, the correlatiou fuuctiou at mid¬ 
rapidity {y = 0 ) is giveu by [ 22 ]: 

CF((j)) — ^ f _ dNgss ( 0 ) \ 

Nmax \PTdpTd(j)dy prdprdi^dy) 

where Nmax uormalizes the correlatiou. Iu this study the value of = 1 is cousidered with lower euergy cutoff 
at 7.6. 

We propose that by measuriug the TPC fuuctiou, oue cau iuvestigate the existeuce of the Mach coue iuter- 
actiou. Heuce, usiug our uumerical data, the TPC fuuctiou for two iuteractiug Mach coues is compared to its 
uou-iuteractiug couuterpart (see Fig. [ 6 ]). The TPC fuuctiou for uou-iuteractiug Mach coues is calculated by 
summiug up the TPC fuuctiou of each Mach coue iu the abseuce of the other oue. Fig. [ 6 ^ shows that prior to 
the iuteractiou, both curves are iu good agreemeut, while after the iuteractiou (see Fig. ( 6 ) 3 ), the TPC fuuctiou 
for the case of iuteractiug Mach coues differs siguificautly from the uou-iuteractiug case. This shows that if 
the Mach coues shock waves do uot iuteract with each other, the TPC fuuctious for the uou-iuteractiug aud 


( 12 ) 
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(a) 



(b) 

Figure 6. Two particle correlation (TPC) function for interacting and non-interacting Mach 
cones for the case T^ed = 0-3 in a domain of 500 x 500 at time, (a) t = 180 (before the 
interaction), and (b) t = 720 fm/c (after the interaction). 

interacting cases should remain the same. However, due to the interaction, the functions deviate from each 
other. 


4. Conclusion 

In this paper we have investigated Mach cone interactions and the resulting RM instability in a relativistic 
fluid. For the numerical simulations, the recently developed LB model for relativistic flows has been extended 
to handle the ideal gas equation of state and external forces. Our results show that the interaction of two Mach 
cone shock waves, and in particular the interaction of a shock front of a Mach cone with the density variations 
generated by the other one, leads to the growth of the RM instability in the direction of the advancing shock 
front. Regarding the thermal behaviour of this phenomenon, we have shown that the average temperature of 
the media decreases because of the instability. 

To single out the effect of the instability on this cooling process, we implemented shock tube RM instability 
simulations, which confirm that the instability causes a decrease in the average temperature. Several simulations 
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have been performed for different initial temperatures and density ratios, which demonstrate that decreasing the 
initial temperature and/or increasing the density ratio, enhances the cooling. This is in line with the analytical 
relation for the linear growth rate of the relativistic RM instability, Eq.([8]), since decreasing the temperature 
and/or increasing the density ratio, enhances the instability. 

We have also shown that the interaction of Mach cones significantly affects the TPC. Therefore, comparing the 
observed TPC with the measured TPC of non-interacting Mach cones, may provide a new tool to experimentally 
identify the interaction. The results of this paper may be relevant to phenomena characterised by the presence 
of Mach cones and/or RM instability in astrophysics, high energy physics, and plasma physics. 

We acknowledge financial support from the European Research Council (ERC) Advanced Grant 319968-FlowCCS. We 
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